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INTRODUCTION 


This  is  the  final  report  on  the  research  contract,  “The  Application  of  Inverse 
Methods  to  the  Ocean  Environment,”  for  the  period,  November  1,  1994  through 
October  31,  1995. 

Under  support  of  this  and  other  research  funding,  we  have  developed  and  imple¬ 
mented  an  inversion  theory  for  imaging  reflectors  in  a  host  medium  and  estimating 
the  changes  in  medium  parameters  across  those  reflectors. 

This  is  a  research  project  in  the  Center  for  Wave  Phenomena,  an  interdisciplinary 
educational  and  research  program  at  the  Colorado  School  of  Mines;  the  principal 
investigator  is  the  director  of  that  center. 

Technology  Transfer 

A  major  source  of  funding  for  the  Center  for  Wave  Phenomena  is  a  consortium 
project  supported  by  thirty-two  companies  in  the  oil  industry:  explorers/producers, 
seismic  data  processors  and  computer  companies.  We  would  like  to  believe  that  the 
success  of  our  program  in  this  this  period  of  dramatic  contraction  in  the  industry  we 
serve  is  a  measure  of  the  quality,  relevance  and  potential  applicability  of  our  research 
in  finding  energy  reserves. 

The  Center  also  has  support  from  the  Ocean  Acoustics  Branch  of  ONR,  the  ACTI 
program,  providing  partnership  between  the  energy  industry  and  the  national  labs 
(Livermore,  Los  Alamos,  Oak  Ridge),  the  Army  Research  Office,  the  Gas  Research 
Institute,  and  the  Foundation  of  the  Society  of  Exploration  Geophysicists. 

Papers 

There  was  one  paper  produced  partially  under  the  support  of  this  project  during 
the  current  year. 

Migration  Velocity  Analysis,  Zhenyue  Liu,  CWP-168,  Ph.  D.  Thesis,  previously  sent 
to  ONR. 

A  brief  excerpt  of  this  paper  was  presented  in  our  consortium  project  review. 
That  excerpt  is  attached  because  it  addresses  the  research  completed  in  the 
current  project  year,  only. 

Test  on  Marmousi  Data  of  Velocity  Analysis  by  Perturbation,  Zhenyue  Liu  and 
Norman  Bleistein,  appearing  in  the  CWP  Annual  Project  Review,  CWP-185. 


Presentations 

Presentations  on  our  research  on  this  project  were  presented  at  the  following 
international  meetings. 
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Velocity  Analysis  by  Perturbation,  by  Liu  and  Bleistein,  Expanded  Abstracts.  Society 
of  Exploration  Geophysicists,  October  1994,  p.  1191.  (Presented  by  Liu.) 

Mathematical  Analysis  of  Residual  Moveout  and  Velocity  Analysis,  Zhenyue  Liu  and 
Norman  Bleistein,  Abstracts,  International  Geophysical  Conference  and  Exposi¬ 
tion,  St.  Petersburg,  Russia,  July  1995.  (Presented  by  Bleistein.) 

Mathematical  Analysis  of  Residual  Moveout  and  Velocity  Analysis,  Zhenyue  Liu 
and  Norman  Bleistein,  Workshop  on  Inverse  Methods,  Institute  for  Mathemat¬ 
ical  Analysis,  University  of  Minnesota,  March,  1995.  (Presented  by  Bleistein.) 

Mathematical  Analysis  of  Residual  Moveout  and  Velocity  Analysis,  Zhenyue  Liu  and 
Norman  Bleistein,  Workshop  on  Inverse  Methods,  International  Meeting  of  the 
Society  of  Exploration  Geophysicists,  Houston,  October,  1995.  (Presented  by 
Bleistein.) 

This  talk  was  also  a  “fill-in”  for  an  absent  speaker  at  a  workshop  on  inversion 
at  the  International  Congress  of  Industrial  and  Applied  Mathematics  in  Hamburg, 
Germany,  July,  1995. 


BACKGROUND 

As  noted  above,  this  research  is  part  of  a  larger  integrated  research  project  on 
the  application  of  inverse  methods  to  map  the  interior  of  the  earth  with  sources 
on  land  or  in  the  ocean.  As  such,  we  concentrate  on  modeling  and  inversion  in 
acoustic  and  elastic  media  (although  our  methods  are  just  as  valid  in  electromagnetic 
media).  The  emphasis  is  on  inversion  of  data  from  multi-source/multi-receiver  arrays 
to  produce  a  reflector  map  of  the  probed  medium  and  to  estimate  the  changes  in 
medium  parameters  across  those  reflectors.  Since  our  concentration  is  on  reflection, 
refraction  and  diffraction,  are  methods  are  naturally  high  frequency  methods. 

Inverse  problems  are  inherently  nonlinear.  The  wave  equation (s)  governing  the 
forward  problem  involve  products  of  the  unknown  medium  parameters  and  the  un¬ 
known  propagating  fields  in  the  interior  of  the  medium.  Hence,  it  is  necessary  to 
solve  these  problems  by  iterative  and/or  recursive  methods.  Our  methods  back  prop¬ 
agate  the  observed  data  into  the  unknown  medium,  using  “background”  propagation 
parameters  in  this  process.  As  new  information  about  the  unknown  parameters  is 
gained,  the  background  can  be  revised  to  produce  a  (presumably)  more  accurate 
picture  of  the  unknown  medium.  Over  the  years,  we  have  extended  our  methods 
through  a  hierarchy  of  progressively  more  realistic  source/receiver  configurations  and 
progressively  more  complex  background  media. 

The  primary  relevance  of  our  work  to  ONR  is  the  application  to  mapping  the 
seabed,  seamine  detection,  and  nondestructive  evaluation — imaging  of  flaws  in  solids — 
where  the  method  is  known  as  synthetic  aperture  focusing  technique  or  SAFT.  In  the 
seismic  inverse  problem,  when  the  sources  and  receivers  are  buried  (located  in  wells), 
the  method  is  called  diffraction  tomography. 
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A  major  stumbling  block  for  all  inversion  techniques  is  getting  good  background 
propagation  speeds  because  the  observed  data  on  the  upper  surface  must  be  “propa¬ 
gated”  back  into  the  subsurface  using  this  background.  Poor  estimates  of  background 
speed  yield  misplaced  reflectors  and  poor  estimates  of  the  change  in  medium  param¬ 
eters  across  them.  Furthermore,  amplitude  control  in  an  absolute  sense  is  extremely 
impractical;  amplitude  control  in  a  relative  sense — from  one  experiment  to  the  next — 
is  a  more  realistic  goal.  Thus,  there  is  a  need  for  a  method  to  estimate  background 
velocities  that  does  not  require  precise  amplitude  information. 

The  interest  in  the  oil  industry  for  this  technique  stems  from  the  need  to  properly 
place  possible  oil  or  gas  prospects.  For  seabed  mapping,  this  type  of  velocity  analysis 
is  crucial  to  proper  mapping  of  the  seabed  in  regions  of  complex  structure.  Sonar 
signals,  for  example,  can  “tunnel”  through  the  shallow  part  of  the  seabed,  arriving 
back  at  the  receiver  location  before  reflections  that  have  propagated  totally  within 
the  ocean.  Thus,  a  knowledge  of  the  seabed  in  regions  of  complex  structure  is  a 
necessary  requirement  for  proper  interpretation  of  such  returns. 

In  flaw  characterization  there  is  a  similar  need  for  proper  position  and  sizing  of 
cracks  and  inclusions  in  solids.  Similarly,  more  accurate  location  and  imaging  of 
tumors  and  other  anomalies  in  human  tissue  would  provide  a  dramatic  improvement 
in  non-invasive  analysis  of  medical  problems. 

Another  important  military  application  is  the  location  of  sea  mines.  While  the 
location  of  floating  mines  might  not  require  the  precision  velocity  analysis  of  our 
research,  the  location  of  mines  that  have  sunk  or  been  buried  in  the  seabed  will 
require  this  kind  of  precision.  While  the  “reflector”  in  this  case  may  be  more  like  a 
point  scatterer,  it  requires  only  a  minor  change  in  our  processing  algorithms  to  image 
point  scatterers  rather  than  reflectors.  Such  modifications  have  been  described  in 
Bleistein  [1988],  partially  supported  by  the  Mathematics  Branch  of  ONR  under  this 
program. 

With  the  graduate  student  partially  supported  by  this  project,  the  principal  in¬ 
vestigator  has  been  working  on  various  approaches  to  “velocity  analysis,”  a  technique 
for  obtaining  a  reasonable  accurate  (low  wave  number)  propagation  speed  to  use  in 
the  computer  implementation  of  the  inversion.  As  noted  above,  we  specifically  use 
methods  that  do  not  rely  on  accurate  amplitude  in  the  data  acquisition.  Thus,  our 
methods  use  only  the  reflector  mapping  property  of  our  inversion  and  not  on  its  pa¬ 
rameter  estimation  capability.  In  the  seismic  exploration  literature,  methods  that 
provide  reflector  maps  are  called  “migration”  rather  than  inversion,  although  the  dis¬ 
tinction  between  the  two  terms  has  blurred  in  recent  years.  However,  in  our  papers, 
we  often  refer  to  migration  rather  than  inversion  to  emphasize  to  the  reader  that  this 
component  of  our  more  general  inversion  techniques  can  be  carried  out  with  standard 
migration  algorithms  that  have  little  regard  for  amplitude  accuracy.  Often,  these 
have  the  advantage  of  requiring  less  CPU,  although  our  inversion  and  comparable 
Kirchhoff  migration  codes  have  similar  run  times. 
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Briefly,  the  idea  behind  our  method  is  as  follows.  It  depends  on  the  fact  that  there 
is  multiple  coverage  of  the  unknown  medium  with  a  suite  of  experiments  of  one  of 
the  following  types: 

1.  common  shot  experiment — fixed  source  and  a  line  of  receivers; 

2.  common  offset  experiment — fixed  distance  between  source  and  receiver,  many 
experiments. 

Note  that  for  the  latter,  one  need  only  reorder  the  data  in  the  suite  of  experiments 
of  the  former  type.  The  former  type  is  more  typical.  A  ship  with  a  towed  array 
of  receivers  sets  off  a  source — dynamite  or  mechanical — at  uniform  intervals  as  it 
moves  in  a  line  in  the  ocean.  The  latter  type  is  more  popular  for  processing  because 
the  suite  of  traces  at  different  midpoints  and  one  offset  illuminates  the  entire  survey 
region  while  the  suite  of  traces  at  one  source  and  all  receivers  would  illuminate  only 
a  small  portion  of  the  medium  below  the  receiver  array. 

Suppose  that  the  background  propagation  speed  is  known  down  to  a  particular 
reflector  and  the  objective  is  to  get  a  “best  guess”  for  the  average  background  in  the 
layer  between  that  level  and  a  next  reflector.  Choose  a  transverse  position  (trace 
location)  where  the  reflector  is  imaged  by  more  than  one  common  shot  inversion  or 
by  more  than  one  common  offset  inversion.  For  each  inversion,  the  illumination  will 
require  different  refractions  at  the  upper  interfaces  and  different  traveltimes  to  the 
reflection  point  under  investigation.  Thus,  the  output  of  that  trace  will  show  the 
reflection  point  in  question  at  the  same  depth  for  all  experiments,  only  if  the  “right” 
propagation  speed  is  used  above  the  reflector;  for  the  wrong  propagation  speed,  the 
different  outputs  will  show  the  “same”  reflection  point  at  different  depths.  It  is 
the  curvature  of  this  ensemble  of  erroneous  depths  that  we  exploit  to  determine  the 
velocity  recursively  through  progressively  deeper  layers  in  the  Earth.  Note  that  this 
technique  depends  upon  a  certain  degree  of  coherency  in  this  common  image  gather 
of  outputs;  it  is  necessary  to  “see”  the  images  of  a  common  reflection  point  and 
distinguish  them  from  the  clutter  of  multiple  images  and  noise.  When  the  velocity 
is  determined  in  a  particular  layer,  the  data  is  processed  through  that  layer  to  the 
next  reflector  and  the  method  is  repeated,  either  below  the  new  reflector  or  laterally, 
until  the  entire  medium  is  covered  by  blocks  in  which  the  velocity  analysis  has  been 
carried  out  and  a  reflector  map  has  been  produced.  The  method  within  a  layer  is 
iterative,  correcting  the  curvature  of  the  multiple  images  of  a  common  reflection  point 
until  they  agree  within  some  error  criterion;  the  method  is  recursive  in  that  it  moves 
progressively  deeper  into  the  host  medium  starting  from  the  upper  surface. 

A  complete  discussion  of  this  iteration/recursion  procedure  and  the  interface 
smoothing  can  be  found  in  Liu  and  Bleistein  [1991]  and  Liu  [1995],  and  is  described 
in  more  detail  in  previous  renewal  proposals.  References  to  related  work  can  also  be 
found  in  those  papers,  as  well. 


5 


More  recently  [Liu  and  Bleistein,  1992],  this  research  project  has  focused  on  the 
analytical  relationship  between  the  apparent  curvature  of  the  erroneous  positions  as 
function  of  the  offset  of  the  source  and  receiver  and  the  error  in  velocity.  For  horizon¬ 
tal  reflectors,  this  technique  does  not  require  the  iterations  of  the  former  method;  one 
reads  off  the  “correct”  velocity  from  the  curvature  of  the  arrival  time  as  a  function  of 
the  offset  between  source  and  receiver.  For  small  lateral  variations  in  velocity,  a  per¬ 
turbation  method  allows  us  to  solve  for  an  approximate  laterally  varying  background 
velocity,  as  well. 

The  analysis  of  [Liu  and  Bleistein,  1992]  has  also  been  extended  to  take  account 
of  velocity  gradients.  This  work  is  described  in  Liu  and  Bleistein,  [1993].  As  part 
of  this  more  recent  analysis,  we  also  provide  sensivity  and  error  estimates  for  this 
type  of  velocity  analysis.  A  shortcoming  of  that  method  was  that  it  still  relied  on  a 
near  hyperbolic  curvature  of  the  ensemble  of  images  of  the  same  depth  point  on  the 
reflector  as  a  function  of  offset. 


Work  in  this  project  year 

The  work  completed  in  the  current  project  year  makes  significant  improvement  on 
the  previous  method  used  to  update  the  velocity  in  the  iterative  part  of  the  algorithm. 
This  is  achieved  by  deriving  an  expression  for  the  functional  derivative  of  image  depths 
with  respect  to  parameters  characterizing  the  mean  velocity  or  mean  components  of 
velocity  gradient  in  a  given  layer  in  the  subsurface.  That  expression  is  described  in 
Liu  and  Bleistein  [1995]  (attached)  and  in  detail  in  the  Ph.D.  thesis,  Liu  [1995]. 

The  most  difficult  test  of  this  data  was  carried  out  on  the  Marmousi  data,  a 
test  data  set  produced  for  the  oil  industry  by  the  Institut  Frangais  du  Petrole  in 
1988  [Versteeg  and  Grau,  1991].  The  physical  model  for  this  data  set  is  shown  in 
Figure  2  in  the  attached  paper.  It  shows  over  200  layers.  This  velocity  analysis 
method  produced  a  background  model  shown  in  Figure  5  containing  only  19  layers. 
This  model  was  produced  in  a  “blind  test,”  assuming  no  a  priori  knowledge  of  the 
model  in  the  test.  Figure  6  shows  a  comparison  of  the  true  velocity  with  the  derived 
velocity  on  three  vertical  lines.  One  can  see  that  the  derived  baickground  velocity 
has  lower  frequency  content  than  the  true  velocity.  This  is  to  be  expected.  Since 
we  are  using  traveltime  information  to  derive  the  background  velocity,  we  are  limited 
to  distinguishing  changes  over  “many”  (at  least  three)  wavelengths  at  the  dominant 
frequency  of  the  data.  However,  one  can  also  see  in  this  comparison  that  the  method 
has  produced  a  velocity  with  a  nonzero  vertical  component  of  gradient  in  some  layers, 
in  part,  as  partial  compensation  for  loss  of  detail.  In  fact,  there  were  cases  of  nonzero 
horizontal  gradient,  as  well. 

The  output  produced  by  inverting  the  Marmousi  data  with  the  background  derived 
by  this  velocity  analysis  is  shown  in  Figure  7,  while  the  inversion  produced  with  the 
true  velocity  model  as  background  is  shown  in  Figure  8. 
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The  only  region  in  which  the  output  of  Figure  8  is  significantly  better  than  that  of 
Figure  7  is  below  the  midpoint  at  6km  at  a  depth  of  2.5km  and  below.  In  this  region, 
this  approach  to  velocity  analysis  has  failed.  The  explanation  for  this  is  provided  in 
Liu  [1995]  and  in  the  attached  paper.  It  is  a  consequence  of  the  fact  that  there  are 
multiple  arrivals  at  this  depth  and  a  lack  of  coherency  in  our  images  of  a  common 
reflector  point  from  data  at  different  offsets.  This  was  first  reported  in  Geoltrain, 
and  Brae  [1991].  This  lack  of  coherency  at  depth  can  be  seen  in  Figure  10.  As  a 
comparison,  Liu  shows  in  Figure  12  that  the  same  lack  of  coherency  occurs  when  one 
uses  the  true  velocity  rather  than  the  derived  velocity,  confirming  that  this  anomaly 
is  not  peculiar  to  our  method,  but  is  inherent  in  the  data.  Thus,  we  need  to  improve 
the  basic  approach  to  obtain  a  better  background  velocity  in  this  one  region. 

The  results  in  the  remainder  of  Figure  7  are  extremely  good,  at  least  as  good  or 
better  than  any  blind  tests  reported  in  a  survey  article  on  processing  of  the  Marmousi 
data  [Versteeg,  R.,  1994]. 

In  the  PROPOSED  RESEARCH  section  of  last  year’s  proposal,  we  set  the  Mar¬ 
mousi  data  set  as  a  benchmark  for  our  inversion.  We  also  noted  last  year  that  we 
needed  a  better  ray  tracing  code  to  achieve  this  end.  That  code  was  developed  under 
other  support.  We  believe  that  we  have  substantially  succeeded  in  meeting  the  goal 
to  develop  a  velocity  analysis  method  to  deal  with  data  of  this  degree  of  complexity, 
although,  as  noted  here,  there  is  still  further  work  to  be  done. 
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